Componentes principales y análisis de imágenes

En este documento se explora la aplicación del análisis de componentes principales sobre imágenes.

Las imágenes han sido obtenidas de A database of facial expressions in younger, middle-aged, and older women and men y solo se han descargado las 72 imágenes que se pueden descargar sin haberse registrado.

Para citar el uso de las imágenes se usa la siguiente referencia

Holland, C. A. C., Ebner, N. C., Lin, T., & Samanez-Larkin, G. R. (2019). Emotion identification across adulthood using the Dynamic FACES database of emotional expressions in younger, middle aged, and older adults. Cognition and Emotion, 33, 245-257. doi:10.1080/02699931.2018.1445981.

Carga de las imágenes

Se utilizarán algunas funciones de la librería imager para la carga y la transformación de las imágenes.

A continuación se carga la librería imager y se indica el directorio donde se encuentran las imágenes:

library(imager)
filenames <- list.files(path = "./img", pattern="*.jpg")

El objeto filenames es una lista que contiene las imágenes. Para cargar una imagen con la funcón load.image() se debe indicar la ruta completa. Para completar el nombre de la ruta se usa la función paste0():

im <- load.image(paste0("./img/",filenames[1]))
plot(im)

Ahora se muestra la imagen en escala de grises con la función grayscale:

plot(grayscale(im))

La dimensión de la imagen es:

dim(im)
## [1] 2835 3543    1    3

Esto quiere decir que el objeto im es una sola imagen de \(2835 \times 3543\) pixeles, con tres canales de color (RGB).

Ahora se utiliza la función resize() para reducir la imagen a diferentes tamaños que corresponden aproximadamente al 1%, 5%, 10% y 20% del tamaño original:

par(mfrow=c(2,2))
plot(resize(grayscale(im),28,36),main="1%")
plot(resize(grayscale(im),143,177),main="5%")
plot(resize(grayscale(im),283,354),main="10%")
plot(resize(grayscale(im),567,709),main="20%")

Leer y transformar todas las imágenes

La librería imager contiene una serie de funciones para leer las imágenes de un directorio y convertirlas luego en un dataframe. Para ilustrar este proceso con un poco más de detalles se usará una función, leer_y_trs_img(), construida para mostrar cada paso de la transformación:

leer_y_trs_img<-function(img_name,path=NULL,x_l=10,y_l=10){
  require(imager)
  img_nombre<-paste0(path,img_name) # completa el nombre de la imagen con la ruta
  imagen<-load.image(img_nombre) # carga la imagen
  img_gris<-grayscale(imagen) # convierte la imagen a escala de grises
  img_escalada<-resize(img_gris,x_l,y_l) # reescala la imagen
  return(img_escalada)
}

A continuación se muestra usa la función para cargar una sola imagen:

x<-leer_y_trs_img(filenames[1],path="./img/",x_l = 57,y_l = 71)
plot(x)

Ahora se procede a cargar todas las imágenes que se indican en filenames usando lapply():

lista_imagenes = lapply(filenames, leer_y_trs_img,path="./img/",x_l = 57,y_l = 71) 

El objeto lista_imagenes es una lista de imágenes. Cada componente es una imagen y se puede utilizar así:

par(mfrow=c(1,2))
plot(lista_imagenes[[2]])
plot(lista_imagenes[[8]])

Ahora se vectorizan las imágenes. Es decir, que vista la imagen como una matriz, sus columnas se ponen una debajo de la otra hasta obtener un vector. Estos vectores luego serán las filas de una matriz de datos donde cada pixel representa una variable y cada imagen representa una observación. El resultado de hacer esto para todas las imágenes es una matriz que tiene tantas filas como imágenes y tantas columnas como pixeles tengan las imágenes. La función as.numeric() aplicada sobre cada imagen devuelve un vector. Se aplica entonces la función as.numeric() sobre cada entrada del objeto lista_imagenes con la función lapply().

imagenes_vectorizadas<-lapply(lista_imagenes, as.numeric)
length(imagenes_vectorizadas[[1]]) # longitud del vector que representa la primera imagen
## [1] 4047
length(imagenes_vectorizadas) # cantidad de imágenes
## [1] 72

El resultado, imagenes_vectorizadas, es una lista de vectores que se concatenan llamando rbind() con do.call():

matriz_imagenes = do.call('rbind', imagenes_vectorizadas)
dim(matriz_imagenes) # dimensión de la matriz resultante
## [1]   72 4047
# as.data.frame.imlist es una alternativa

Estadísticos descriptivos

Al tener una matriz de datos se pueden hacer todo lo que se hace con una matriz de datos. Sin embargo, el gran número de variables hace que utilizar la función summary() o la función pairs() sea poco práctico. Una aproximación puede ser generar estadísticas de resumen para cada pixel y luego ver estas estadísticas en forma de imágenes.

Veamos por ejemplo la imagen media:

imagen_media_vec<-apply(matriz_imagenes,2,mean) # calcula el valor promedio para cada pixel (columna)
length(imagen_media_vec) # longitud de la imagen promedio como vector
## [1] 4047

Ahora convirtamos la imagen promedio vectorizada en una imagen usando la función as.cimg():

imagen_media<-as.cimg(array(imagen_media_vec,dim=c(57,71)))
plot(imagen_media)

De igual manera se puede calcular la imagen que representa las desviaciones estándar para cada pixel:

imagen_sd_vec<-apply(matriz_imagenes,2,sd)
imagen_sd<-as.cimg(array(imagen_sd_vec,dim=c(57,71)))
plot(imagen_sd)

Descartando información poco relevante

La imagen de desviaciones estándares nos muestra varias zonas negras, donde el valor es cero. Estas zonas de baja variabilidad corresponden a pixeles poco informativos en esta muestra de imágenes. El fondo de la imagen, común para todos los individuos es, por ejemplo, una zona de baja variabilidad. Usemos el histograma de la imagen de desviaciones estándares para detectar estas zonas de baja variabilidad:

hist(imagen_sd)

Consideraremos que pixeles con una variabilidad inferior a 0.05 no serán relevantes en los siguientes análisis. La función threshold() asigna un 0 a cada pixel cuya intensidad está por debajo de cierto valor y 1 al resto de pixeles y nos devuelve una imagen. A continuación se usa esta función especificando un valor de umbral de 0.05:

imagen_sd_tr<-threshold(imagen_sd,thr = 0.05)
plot(imagen_sd_tr)

La imagen imagen_sd_tr es una máscara. Las zonas en blanco representan los pixeles que se incluirán en los análisis y las zonas en negro los pixeles que se descartarán por su baja variabilidad.

Veamos cuantos pixeles se descartan:

table(imagen_sd_tr)
## imagen_sd_tr
## FALSE  TRUE 
##   739  3308

Calculemos el porcentaje de información descartada:

table(imagen_sd_tr)[1]/sum(table(imagen_sd_tr))*100
##    FALSE 
## 18.26044

Esto quiere decir que aproximadamente el 18% de los datos están en zonas de baja variabilidad. Ahora identificamos los pixeles por fuera de las zonas de baja variabilidad usando la función which() y almacenamos la ubicación de estos pixeles en el objeto mascara:

mascara<-which(imagen_sd_tr==1)

Ahora se extraen los datos correspondientes a los pixeles por fuera de la zona de baja variabilidad usando el objeto mascara para identificar las columnas relevantes en la matriz_imagenes. El resultado se almacena en el dataframe datos_con_mascara:

datos_con_mascara<-matriz_imagenes[,mascara]
dim(datos_con_mascara)
## [1]   72 3308

Análisis de componentes principales

A continuación se procederá a hacer una análisis de componentes principales sobre las imágenes. El primer paso será centrar y escalar la matriz de datos con la función scale(). Esto quiere decir que a cada columna se le resta su media y se divide por su desviación estándar. El resultado se almacena en el objeto datos_mascara_centrados:

datos_mascara_centrados<-scale(datos_con_mascara,center = TRUE,scale = TRUE)
dim(datos_mascara_centrados)
## [1]   72 3308

Notemos que nuestra matriz de datos tiene muchas más variables (3308) que observaciones (72). Esto hace que no se pueda usar la función princomp() para hacer el análisis de componentes principales. En su lugar se usa la función prcomp(), que utiliza la descomposición en valores singulares. A lo sumo se tendrán tantas componentes como observaciones (o filas de la matriz de datos):

modelo_pca<-prcomp(datos_mascara_centrados)
summary(modelo_pca)
## Importance of components:
##                            PC1     PC2      PC3      PC4      PC5     PC6
## Standard deviation     26.5193 25.1164 16.40629 15.33239 12.44363 9.93699
## Proportion of Variance  0.2126  0.1907  0.08137  0.07106  0.04681 0.02985
## Cumulative Proportion   0.2126  0.4033  0.48467  0.55573  0.60254 0.63239
##                            PC7     PC8    PC9    PC10   PC11    PC12   PC13
## Standard deviation     8.96227 8.51391 8.3537 7.88131 7.3197 6.72317 6.4567
## Proportion of Variance 0.02428 0.02191 0.0211 0.01878 0.0162 0.01366 0.0126
## Cumulative Proportion  0.65667 0.67858 0.6997 0.71846 0.7347 0.74832 0.7609
##                           PC14    PC15   PC16    PC17    PC18    PC19    PC20
## Standard deviation     6.11721 5.98341 5.6654 5.55175 5.29583 5.22120 4.97269
## Proportion of Variance 0.01131 0.01082 0.0097 0.00932 0.00848 0.00824 0.00748
## Cumulative Proportion  0.77223 0.78305 0.7928 0.80207 0.81055 0.81879 0.82627
##                           PC21    PC22    PC23    PC24    PC25    PC26    PC27
## Standard deviation     4.87272 4.76576 4.69384 4.58053 4.47333 4.32228 4.26341
## Proportion of Variance 0.00718 0.00687 0.00666 0.00634 0.00605 0.00565 0.00549
## Cumulative Proportion  0.83345 0.84031 0.84697 0.85331 0.85936 0.86501 0.87051
##                           PC28    PC29    PC30    PC31    PC32    PC33    PC34
## Standard deviation     4.20166 4.13948 4.09074 4.01823 3.93728 3.91091 3.74792
## Proportion of Variance 0.00534 0.00518 0.00506 0.00488 0.00469 0.00462 0.00425
## Cumulative Proportion  0.87584 0.88102 0.88608 0.89096 0.89565 0.90027 0.90452
##                           PC35    PC36    PC37    PC38    PC39    PC40   PC41
## Standard deviation     3.71081 3.67334 3.66160 3.55966 3.47681 3.46320 3.3525
## Proportion of Variance 0.00416 0.00408 0.00405 0.00383 0.00365 0.00363 0.0034
## Cumulative Proportion  0.90868 0.91276 0.91681 0.92064 0.92430 0.92792 0.9313
##                           PC42    PC43    PC44    PC45    PC46    PC47    PC48
## Standard deviation     3.32706 3.29960 3.18259 3.14377 3.10711 3.06952 3.06146
## Proportion of Variance 0.00335 0.00329 0.00306 0.00299 0.00292 0.00285 0.00283
## Cumulative Proportion  0.93467 0.93796 0.94102 0.94401 0.94693 0.94978 0.95261
##                           PC49    PC50   PC51    PC52    PC53    PC54    PC55
## Standard deviation     3.00845 2.95087 2.9350 2.89881 2.88403 2.84588 2.82912
## Proportion of Variance 0.00274 0.00263 0.0026 0.00254 0.00251 0.00245 0.00242
## Cumulative Proportion  0.95534 0.95798 0.9606 0.96312 0.96564 0.96808 0.97050
##                           PC56    PC57    PC58    PC59    PC60    PC61    PC62
## Standard deviation     2.76507 2.74264 2.69070 2.65895 2.60921 2.60380 2.52446
## Proportion of Variance 0.00231 0.00227 0.00219 0.00214 0.00206 0.00205 0.00193
## Cumulative Proportion  0.97281 0.97509 0.97728 0.97941 0.98147 0.98352 0.98545
##                           PC63    PC64    PC65    PC66    PC67    PC68    PC69
## Standard deviation     2.51638 2.48671 2.43214 2.38532 2.33871 2.28619 2.24074
## Proportion of Variance 0.00191 0.00187 0.00179 0.00172 0.00165 0.00158 0.00152
## Cumulative Proportion  0.98736 0.98923 0.99102 0.99274 0.99439 0.99597 0.99749
##                           PC70    PC71      PC72
## Standard deviation     2.14099 1.92727 1.178e-14
## Proportion of Variance 0.00139 0.00112 0.000e+00
## Cumulative Proportion  0.99888 1.00000 1.000e+00

Las primeras 33 componentes principales explican aproximadamente el 90% de la variabilidad de los datos. Tomaremos los primeros 33 vectores propios (o en este caso vectores singulares) para formar las 33 primeras componentes principales:

eigen_modes_prin<-modelo_pca$rotation[,1:33]
dim(eigen_modes_prin)
## [1] 3308   33

Ahora se multiplicará cada imagen por los 33 vectores propios para representarla usando solo 33 variables:

datos_red<-datos_mascara_centrados%*%eigen_modes_prin
dim(datos_red)
## [1] 72 33

El resultado es una matriz de 72 filas y 33 columnas. Cada fila corresponde a un individuo y cada columna a una de las 33 primeras componentes principales.

Análisis de características en el espacio de componentes principales

Trataremos de ver si las en el espacio de las componentes principales es posible observar agrupaciones de imágenes por características. En la base de imágenes los individuos pueden ser jóvenes, adultos o adultos mayores, hombres o mujeres o tener una expresión facial específica: neutral (n), disgusto (d), miedo (f), felicidad (h), tristeza (s) o rabia (a).

Esta información esta codificada en el nombre de la imagen, veamos cómo:

head(filenames)
## [1] "004_o_m_a_a.jpg" "004_o_m_a_b.jpg" "004_o_m_d_a.jpg" "004_o_m_d_b.jpg"
## [5] "004_o_m_f_a.jpg" "004_o_m_f_b.jpg"

Así, las características pueden extraerse usando la función substr():

edad<-substr(filenames,start=5,stop=5)
sexo<-substr(filenames,start=7,stop=7)
expresion<-substr(filenames,start=9,stop=9)

Veamos ahora cómo se agrupan las características en el espacio de las primeras cinco componentes principales usando la función pairs():

pairs(datos_red[,1:5],col=as.factor(edad), main="Edad")

pairs(datos_red[,1:5],col=as.factor(sexo), main="Sexo")

pairs(datos_red[,1:5],col=as.factor(expresion), main="Expresión")

¿Qué significa cada modo?

Veamos ahora la representación de cada vector propio. Para ello convertiremos cada vector propio en una imagen. Primero se crea una matriz que tendrá 33 filas, una para cada vector propio y \(57\times 71\) columnas. Usando la máscara definida anteriormente se llenan los pixeles (columnas) que sí tienen variabilidad:

modas<-matrix(0,ncol=57*71,nrow=33)
modas[,mascara]<-t(eigen_modes_prin)

La función sweep() permite reescalar las modas multiplicándolas por las desviaciones estándar y luego sumándole la imagen media:

modas_escaladas<-sweep(modas,2,FUN='*',imagen_sd_vec)
modas_centradas<-sweep(modas_escaladas,2,FUN='+',imagen_media_vec)

Ahora se convierte cada moda en una imagen:

modas_img_escalanat<-list()
modas_img<-list()
for (i in 1:33){
  modas_img[[i]]<-as.cimg(modas[i,],x=57,y=71)
  modas_img_escalanat[[i]]<-as.cimg(modas_centradas[i,],x=57,y=71)
}

Finalmente visualizamos las primeras seis modas (sin reescalar):

par(mfrow=c(3,2))
lapply(X=1:6,FUN=function(x){plot(modas_img[[x]],main=paste0("Moda ",x))})

## [[1]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1 
## 
## [[2]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1 
## 
## [[3]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1 
## 
## [[4]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1 
## 
## [[5]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1 
## 
## [[6]]
## Image. Width: 57 pix Height: 71 pix Depth: 1 Colour channels: 1

Actividad

Construya un clasificador para las expresiones faciales usando regresión logística y análisis de componentes principales.

FIN